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Organization of the magnetosphere during substorms 

^1 1 

T. Zivkovic, and K. Rypdal, 

Abstract. The change in degree of organization of the magnetosphere during substorms 
is investigated by analyzing various geomagnetic indices, as well as interplanetary mag- 
netic field a^component and solar wind flow speed. We conclude that the magnetosphere 
self-organizes globally during substorms. but neither the magnetosphere nor the solar wind 
become more predictable in the course of a substorm. This conclusion is based on anal- 
ysis of five hundred substorms in the period from 2000 to 2002. A minimal dynamic-stochastic 
model of the driven magnetosphere that reproduces many statistical features of substorm 
indices is discussed. 



1. Introduction 

The complexity of the magnetosphere has been exten- 
sively studied and one of its descriptions is based on the 
paradigm of self-organized criticality (SOC) coined by Bak 
et al. [1987]. Systems that exhibit SOC activate many de- 
grees of freedom, their interactions are local, but due to the 
excitation of a scalo-frcc hierarchy of avalanches when the 
system is slowly driven to a criticality threshold, long-range 
interactions develop, and the dynamics on different spatial 
and temporal scales is essentially the same. It has been 
shown that the magnctosphcric-ionospheric system exhibit 
signatures characteristic of SOC-dynamics [Voros, 1991; 
Uritsky and Pudovkin, 1998], and models for such descrip- 
tion have been developed [Valdivia et al, 2006; Klimas et 
al., 2000; Chapman et al., 1998; Chang, 1999; Uritsky et 
al., 2002; Kozelov and Kozelova, 2003]. On the other hand, 
the magnetospheric-ionospheric system has also some sig- 
natures of intermittent turbulence [Golovchanskaya et al., 
2008], non-equilibrium phase-transitions [Sitnov et al., 2001] 
as well as low-dimensional chaos [Sharma ct al., 1993]. It has 
been recently reported in Zivkovic and Rypdal [2011] that 
the magnetosphere, interplanetary magnetic field (IMF) z- 
component Bz and solar wind flow speed v become more pre- 
dictable and more persistent during magnetic storms, while 
only the magnetosphere reduces the effective number of de- 
grees of freedom through self-organization. In that study 
the magnetosphere was studied through analysis of the ge- 
omagnetic indices SYM-H and Dst, since it is known that 
these indices respond to the intensification of the ring cur- 
rent during magnetic storms [Wanliss, 2006]. In this article 
we analyze how different magnctosplicric indices respond to 
the much more frequent and short-living events called sub- 
storms. 

Different geomagnetic indices represent different parts of 
the magnetosphere and respond to different dynamics. Ge- 
omagnetic indices studied in this article, are downloaded 
from World Data Center, with 1-min resolution. The most 
commonly used index for substorm studies is the auroral 
electrojct index (AE) defined as the difference between the 
AU index, which measures the strength of the eastward elec- 
trojct in the auroral zone, and the AL index, measuring the 
westward electrojct current, and is usually derived from 12 
magnetometers positioned below the auroral oval [Davies 
and Sugiura, 1966]. We also use minute data for the IMF 
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component Bz, as well as minute data for the solar wind 
bulk velocity v along the Sun-Earth axis. They are both 
retrieved from the OMNI satellite database and are given in 
GSE coordinate system. 

We have not analyzed other plasma parameters, since the 
plasma instruments can be influenced by solar X rays or en- 
ergetic particle precipitation and are more unstable than 
the magnetometers. We also analyze the polar cap mag- 
netic activity (PC), as well as the AU and the AL index. 
During magnetically disturbed times the westward electro- 
jct, whose proxy is the AL index, increases abruptly due to 
currents from the magnctotail. On the other hand, the east- 
ward electrojct, whose proxy is the AU index, increases due 
to the partial ring current closure via the ionosphere in the 
evening sector [Feldstein et al, 2006]. Therefore, through 
the analysis of the AL and AU index we can get insight into 
the dynamics of different parts of the magnetosphere. The 
PC index monitors geomagnetic activity over the polar caps 
caused by changes in the IMF and the solar wind. This index 
is mostly influenced by fleld-aligned currents which flow at 
the poleward rim of the auroral oval, and is also sensitive to 
the ionospheric Hall currents in the polar cap [Vennestr0m 
et al, 1991], which are particularly dominant in the summer 
time. Since field-aligned currents are closely related to the 
auroral electrojets, the linear correlation of the PC and AE 
indices is of the order of 0.8-0.9 [Vennestr0m et al., 1991]. 
Here we use the northern polar cap index measured at the 
Danish geomagnetic observatory in Thule (86.5°N). 

Magnetic substorms are associated with release and stor- 
age of energy and momentum from the solar wind to the 
magnetosphere. They consist of three phases. In the growth 
phase, typically lasting for about one hour, loading of the 
magnetic flux and energy into the magnctotail takes place. 
It is succeeded by the expansion phase (substorm onset), 
lasting 30-60 minutes, when hot plasma "unloads" earth- 
ward, leading to sudden brightenings of the polar aurora, 
and plasmoids are ejected away from Earth in the far tail. In 
this phase dissipation is dominant and formation of a storm 
ring current can take place. In the ionosphere substorm 
onset is characterized by a rise of the westward electrojct 
current, which forms the substorm current surge with field 
aligned currents. The recovery phase returns the magneto- 
sphere to its quiet state. The duration of this this phase is 
1-2 hours. 

Data about the times of substorm onsets arc found in 
Frey et al. [2004], whose database is from the period be- 
tween 2000 and 2002. These substorms were detected by 
the FUV instrument on the IMAGE spacecraft. Observa- 
tions covered the peak of the last solar cycle. According 
to Frey et al. [2004] a substorm onset is only accepted as a 
separate event if at least 30 minutes have passed after the 
previous event. It is also required that local brightening of 
aurora occurs and that the aurora expands to the poleward 
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boundary of the auroral oval and spreads azirriuthally in lo- 
cal time for at least 20 minutes. The latter criterion excludes 
pseudo-breakups which do not develop into full substorms. 

The remainder of the paper is organized as follows: sec- 
tion 2 gives a brief overview of the methods used in the 
analysis of our data. Section 3 shows the results from appli- 
cation of these methods to the substorm data and section 4 
is reserved for discussion. 

2. Methods 

2.1. Recurrence-plot analysis 

Recurrence-plot analysis was developed by Eckmann et 
al, [1987] and is very useful in studies of short and non- 
stationary time scries. A comprehensive review of the 
method and its applications can be found in Marwan et al. 
[2007]. Substorm durations are at most a few hours and 
all indices show non-stationary behavior during the events. 
Recurrence-plot analysis is very suitable for handling such 
short non-stationaxy time series. The method is useful for 
low-dimensional deterministic dynamical systems as well as 
for high-dimensional systems and for stochastic signals. It 
can also provide useful information about non-autonomous 
systems. An impression of the versatility of this technique 
can be obtained from the special issue edited by Marwan et 
al. [2008] . For low-dimensional systems 

^ = fm,t], (1) 

recurrence plots are based on the recurrences of a trajectory 
z(t) on the d-dimcnsional attractor in a p > rf-dimensional 
phase space. If the system is autonomous, i.e. no explicit 
time dependence of the phase-space flow f[z], and if the at- 
tractor of trajectory has dimension d, Takens' time-delay 
method [Takens, 1981] can be used to construct an m > 2d- 
dimensional embedding space on which the attractor can be 
mapped continuously and one-to-one. The embedding space 
is constructed from the time series Xi = g{z[i5t]), where 
g{z) is the measurement function. Here t = i6t, and 5t is 
the sampling time of the time series. The mapping is given 
by 

Xj = (iCj, Xi-\-T, . - - , iCj-|-(7^_i)T-)? (2) 

where r is the time delay. There are practical constraints on 
useful choices of the time delay r. If r is much smaller than 
the autocorrelation time the image of the attractor in the 
embedding space becomes essentially one-dimensional. If r 
is much larger than the autocorrelation time, noise may de- 
stroy the deterministic connection between the components 
of x(f), such that our assumption that z(t) determines x(t) 
will fail in practice. A common choice of r has been the 
first minimum of the autocorrelation function, but it has 
been shown that better results are achieved by selecting the 
time delay as the first minimum in the average mutual in- 
formation function [Abarbanel, 1996], which we also use in 
this article. 

The recurrence-plot analysis deals with the trajectories 
in the embedding space. If the original time series x{t) has 
N elements and a time delay r, we have a time series of 
N-{m- 1)t vectors x(t) for t = 1, . . . , - (m - 1)t. This 
time series constitutes the trajectory in the reconstructed 
embedding space. 

The next step is to construct a N—{m—l)TxN—{m—l)T 
matrix Rij consisting of elements and 1. The matrix el- 
ement is 1 if the distance is )|xj — Xj)| < e in the 
reconstructed space, and otherwise it is 0. The recurrence 
plot is simply a plot where the points for which the 
corresponding matrix element is 1 is marked by a dot. For 
a deterministic system the radius e is typically chosen as 



small fraction of the diameter of the reconstructed attrac- 
tor, but varies for different sets of data. In our analysis we 
have used 10% of the total extent of the set spanned by 
the trajectory analyzed in the embedding space. This is a 
rule of the thumb generally accepted in the recurrence-plot 
literature [Marwan et al., 2007]. The results are not very 
sensitive to this choice, but issues arise if e is too large (to 
crude resolution and practically no information) or too small 
(poor statistics). 

Dynamical systems with a large number of independent or 
weakly dependent degrees of freedom can only be described 
either by large-scale numerical simulation or by stochastic 
methods. For such systems the phase-space attractor is also 
high-dimensional and cannot be mapped one-to-one onto a 
low-dimensional time-delay embedding space. Nevertheless, 
the evolution of the "projection" of the phase-space vector 
onto the embedding space usually provides valuable informa- 
tion if the measurement function is carefully chosen to be 
sensitive to variation of those degrees of freedom that are 
in focus of our interest (for instance we should use the AE 
index rather than the Dst index as measurement function if 
substorm activity is studied). Some of this information can 
be discerned from the recurrence plots, even if the recurring 
states are not recurrences of the full phase-space vector, but 
only of the projections. This is true also if the full set of dif- 
ferential equations is non-autonomous. The time-dependent 
forcing invalidates the embedding theorem, but we arc not 
assuming a perfect embedding anyway. We just have to be 
aware that dynamical features we observe in the recurrence 
plot may well reflect the dynamics of the forcing and not 
only the internal dynamics of the unforced system. 

Of particular interest is analysis based on the diagonal 
line structures of the recurrence plots. If the projected tra- 
jectory has a tendency to repeat its path when it recurs to 
the same region in embedding space, recurrences tend to 
be localized along unbroken diagonal segments. The longer 
these segments are, the more predictable is the path. We 
define the average inverse diagonal line segment length F: 

r^r^) = ^r^p(z)/^p(0, (3) 
I I 

where P{V) is a histogram over diagonal lengths: 

N l-l 

P{l) = ^ (1 - -Ri-i,i-i)(i - Ri+i,3+i) n Ri+k,i+k.{^) 

i,j = l k=0 

It was shown heuristically in Zivkovic and Rypdal [2011] that 
r can be used as a proxy for the largest Lyapunov exponent 
in deterministic systems and as a measure of persistence in 
stochastic systems. For a deterministic system we demon- 
strated the connection with the largest Lyapunov exponent 
A by computing A and F along a period-doubling route to 
chaos for the Lorcnz system. By computing F and the self- 
similarity exponent h for rmmerically generated fractional 
Brownian motions (fBm) with h in the range < h < 1 
we established the relation F = 0.72 — 0.57/i in the range 
of persistent Brownian motions (0.5 < h < 1). An fBm 
is a Gaussian stochastic process {x{t)} which satisfies the 

self-similarity condition x{t) = \^'^x{\t) for all A. Here = 
denotes identity in distribution. For a given time-series h 
can be estimated from the variogram [Zivkovic and Rypdal, 
2011]. For a chaotic, deterministic system the inverse of the 
largest Lyapunov exponent is a measure of predictability. 
For a persistent stochastic process h is also such a measure, 
and F decreases with increasing h. Thus we suggest use the 
inverse of F as a measure of predictability. Since signals of 
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interest are a mixture of deterministic and stochastic compo- 
nents, and many stochastic processes are neither self-similar 
nor Gaussian, we do not consider F as an absolute measure 
of predictability, but when a decrease of F takes place in a 
time series under certain conditions, e.g. during a magnetic 
storm, we interpret this as increase in predictability of the 
dynamics under these conditions. 

We compute F for embedding dimension m = 1. An 
obvious advantage of this choice is that we don't have to 
worry about the choice of time delay r. A higher value of 
m does not make much more sense since the dynamics of 
the magnetosphere and the solar wind is strongly influenced 
by a high-dimensional/stochastic component and cannot be 
unfolded in any higher-dimensional embedding space. Also, 
if the embedding dimension is inappropriately high, spuri- 
ous long diagonal lines will appear in the recurrence plot. 
Moreover, the way F depends on other predictability mea- 
sures like Lyapunov exponent or self-similarity exponent is 
rather insensitive to the choice of m, hence F can be used 
to detect changes in predictability irrespective of the choice 
of embedding dimension. 

2.2. A test for determinism 

We shall adopt a terminology where a physical system is 
deterministic if it can be described as a low-dimensional dy- 
namical system, i.e. by the autonomous version of equation 
(1). A test of determinism was developed by Kaplan and 
Glass [1992, 1993] which takes advantage of the fact that 
the trajectory through a given position in phase space is 
completely determined by this position, and if a trajectory 
recurs to the vicinity of this point the tangents of the tra- 
jectory in these two points are approximately parallel. This 
is in contrast to a stochastic system, where the directions 
of the tangents are independent for the two points. Let us 
assume that the phase space is divided into boxes. In our 
test, a box size is chosen as the average distance a phase- 
space point moves in the m-dimensional embedding space 
during one time step. The displacement of the trajectory 
inside a box, in m-dimensional phase space is given from 
the time-delay embedding reconstruction: 



Ax(t) = [x{t + b) - 
x{t + (m ^ 



x{t),x{t + T + b) - x{t + t), . . . , 
^l)r + b)~x{t+{m-l)T)], (5) 



where b is the characteristic time the trajectory spends 
inside a box. In our test b = 1 time step. The tan- 
gent for the kth pass through box j is the unit vector 
Ukj = Axfe,j(t)/|Axfc,j(f)|. The averaged tangent in the 
box is 



Hi ^ — ' 



(6) 



where rij is the number of passes of the trajectory through 
box j. In the case of deterministic dynamics and finite box 
size, Vj will not depend very much on the number of passes 
Uj, and Vj will converge to 1. In contrast, for the trajec- 
tory of a random process with independent increments, Vj 
will decrease with Uj as Vj ^ n- . Thus, to obtain better 
statistics for the description of average tangents we compute 
the average Vj as a function of the number of passes through 
a box: 



of projected trajectories starting at almost the same posi- 
tion in the reduced embedding space. In those cases when 
the system can be described as a low-dimensional dynam- 
ical system (equation (1) without the explicit time depen- 
dence in the flow field), we have I/„ = 1 for all n, provided 
the embedding dimension m is sufficiently large to unfold 
the attractor. In this case F will reflect the magnitude of 
the largest Lyapunov exponent, and the two tests clearly 
measure different properties. For a stochastic process L„ is 
independent of m [Zivkovic and Rypdal, 2011] and indepen- 
dent of its persistence, which is easily demonstrated by com- 
puting it for fBms with varying self-similarity exponents h. 
However, for the systems considered in this paper the signals 
contain both a deterministic and a stochastic component. In 
this case the phase-space trajectory cannot be mapped one- 
to-one onto the embedding space for any choice of m, and 
the distinction between the two methods is less clear. One 
could envisage that an increase of the low-dimensional (de- 
terministic) component relative to the stochastic one could 
enhance both L„ and F~^. On the other hand an increase 
in predictability in either the deterministic or the stochas- 
tic component, without change in relative strength of the 
two components, could also influence both measures. Thus, 
in those cases where both measures either increase or de- 
crease together, it is difficult to decide whether the cause is 
a change in determinism or in predictability. However, in sit- 
uations when only one of the measures undergoes a change, 
or they change in opposite direction, the interpretation is 
unambiguous. 

Kaplan and Glass [1992, 1993] also suggest another test 
of determinism, which does not suffer from this ambigu- 
ity. It takes advantage of how the measure of determin- 
ism, the curves L„, changes when we construct a new sur- 
rogate time series for which the effect of nonlinearity in a 
low-dimensional system has been corrupted. This surrogate 
signal has the same power spectral density (and hence same 
auto-covariance) as the original signal, but the phases of 
the Fourier coefficients have been randomized. From ap- 
plication of this procedure to synthetic stochastic and low- 
dimensional signals we have gained support for the con- 
jecture that randomization of phases do not change these 
curves for neither stochastic or high-dimensional, nor low- 
dimensional, linear systems. Only for low-dimensional, non- 
linear systems will the surrogate data curves lie below those 
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L„ = {V,)„^.=„, (7) 

where this average is done over all boxes with same number 
n of trajectory passes. 

There are obvious similarities between this test for de- 
terminism and the one for predictability developed in the 
previous subsection. Both measure the degree of divergence 



Figure 1. Mean vector length L„ vs. number of pas- 
sages n: bullets are derived from numerical solutions to 
the MG equation, stars are derived from these solutions 
after randomization of phases of Fourier coefficients. In- 
set in the figure shows same characteristics for the fO-U 
process. The error bars are of the same order as the 
symbol size. 
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based on the original data. These features were demon- 
strated for numerical solutions of the Lorenz system (low- 
dimensional, nonlinear and chaotic) in Zivkovic and Rypdal 
[2011]. Here we demonstrate the same for L„ versus n for 
the Mackey-Glass (MG) equation Mackey and Glass [1977]: 



dx , , x(t — 5) 

— = -bxU) + a- ^ 

dt ^ ' l + x{t 



5)'= 



(8) 



Unlike the Lorenz system, the attractor of this system is in 
general high-dimensional. However, for some set of param- 
eters, e.g. a = 0.2,6 = 0.1,6 = 100, c = 10, the attractor 
of MG dynamics can be low-dimensional and L„ derived 
from the MG equation will fall more slowly with increasing 
n than that for the randomized version, as shown in Figure 
1. The inset in the same figure shows the same results for 
a fractional Ornstein-Uhlenbeck (fO-U) stochastic process. 
The stochastic equation used to generate such a process was 
described in [Zivkovtc and Rypdal, 2011]. The coefficients in 
this equation used to generate Figure 1 are generated from 
fitting the variogram of the numerically generated process 
to the variogram of the AE index by means of least square 
regression (for definition of the variogram see equation (12) 
in Zivkovic and Rypdal [2011]). In both fO-U and MG equa- 
tion, the embedding dimension used is m = 8. Realizations 
of a stochastic process like fO-U is indistinguishable from re- 
alizations of a measurement function of a high-dimensional 
deterministic system for which the embedding dimension is 
too small to unfold the attractor. The L„-curve for the 
fO-U process is unchanged after randomization of phases 
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Figure 2. Recurrence plots for tAE: a) for the entire 
day of September 8th, 2002, containing two distinct sub- 
storms whose onsets are marked by arrows, b) Blow-up 
of recurrence plot for second substorm (marked by square 
in panel (a)). 



and shows that this test for determinism is negative for 
stochastic (and high-dimensional) systems. We have veri- 
fied that this is the case also for strongly persistent (highly 
predictable) fO-U processes. This is completely reasonable 
on theoretical grounds; the long-range persistence in an fBm 
only depends on the predominance of low frequencies in the 
power spectrum, not on correlation between the phases of 
the Fourier coefficients. Thus, this test is not measuring 
predictability and is a useful test for detecting changes in 
determinism in time series. 

3. Results 

In Rypdal and Rypdal [2010a] it was shown that the fluc- 
tuation amplitude (or more precisely; the one-timestep in- 
crement) Ay{t) of the AE index is on the average pro- 
portional to the instantaneous value y{t) of the index. 
This gives rise to a special kind of intermittency asso- 
ciated with multiplicative noises, and leads to a non- 
stationary time series of increments. However, the time se- 
ries Ay{t)/y{t) is stationary, implying that the stochastic 
process x{t) = log y{t) has stationary increments. Thus, a 
signal with stationary increments, which still can exhibit a 
multifractal intermittency, can be constructed by consider- 
ing the logarithm of the AE index. We use transformed 
versions of geomagnetic indices: tAE= \og{AE), tAL= 
log(0.1928 |AL| + 10.50), tAU=log(0.1167 |AU[ + 13.50), 
tPC= log(0.0503 |PC| + 0.1978), while B, and v have in- 
crements which are not strongly dependent on their magni- 
tude, and do not need transformation to obtain stationary 
increments. 

It is well known that the AE index exhibits scale-free 
characteristics which often is associated with stochastic dy- 
namics since its power spectral density has two distinct 
power-law regimes [Tsurutani et al, 1990]. However, it has 
also been shown that apart from colored noise which dom- 
inates the dynamics of the AE index on time scales up to 
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Figure 3. a) Ensemble mean of F averaged over 500 
substorms. b) Ensemble mean of 2h. Triangles represent 
tAE, squares tAU, stars tPC, and circles tAL. The stan- 
dard errors for all points are in the range 0.02 < a < 0.04. 
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100 minutes, there are also signatures of low-dimensional, 
chaotic characteristics as well [Athanasiu and Pavlos, 2001]. 
These properties will be discerned from the analysis of de- 
terminism of geomagnetic indices in section 3.2. 

3.1. Predictability analysis 

First, we discuss two magnetic substorms that occurred 
on September 8th, 2002. The first substorm onset was reg- 
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Figure 4. a) Ensemble mean for F. b) Ensemble mean 
for 2h. Stars represent Bz and diamonds v. Standard 
error for F is of the order a « 0.02, while for 2h it is 
a « 0.04. 
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istered at 5:01 UT and the other at 21:26 UT. The second 
substorm onset started after 13 hours of northward IMF B^. 
In Figure 2a the transition between dark and white bands in 
the recurrence plot indicates changes in the dynamics in the 
tAE. The first transition marked by an arrow occurs around 
the first substorm onset. The next appears at the onset of 
the second substorm. For instance, the broad white band in 
Figure 2a corresponding to i > 300, j < 300 indicates that a 
region in embedding space visited before the substorm on- 
set is never visited after onset. Figure 2b contains a blow-up 
of the recurrence plot for the second substorm. The image 
pales away from the main diagonal, which indicates that 
the process is non-stationary due to the intensification and 
the movement of ionospheric currents during the substorm. 
Since the plot pattern changes during the substorm we might 
expect that F, which is defined from the diagonal lines of the 
plot, will change as well. We obtain F for tAE, tAU, tPC, 
tAL, Bz, and v for 500 substorms whose onsets were taken 
from the database of Frey et al. [2004]. The variation of 
F over a 6 hour time interval, three hours before and three 
hours after the substorm onset, is computed from recurrence 
plots derived from thirty-minute windows, such that 12 F- 
values are obtained for each substorm. The time evolution 
of F is computed for all 500 substorms and averaged. The 
results are presented in Figure 3a. Apparently, there is no 
significant variation of the ensemble average of F over the 
duration of a substorm for any of the geomagnetic indices. 

On time scales less than 100 minutes all quantities ana- 
lyzed can be modeled as non-stationary multifractal motions 
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Figure 5. Mean vector length Ln vs. number of pas- 
sages n: bullets are derived from the entire tAE time 
series, stars from this time series after randomization of 
phases of the Fourier coefficiens. Standard errors are of 
the same order as the symbol size. 



Figure 6. L„ vs. n, a) squares are derived from tAE 
during substorms, triangles from the entire tAE time se- 
ries, b) squares from tAU during substorms, triangles 
from entire tAU time series. Standard errors are of the 
same order as the symbol size. 
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whose persistence can be characterized by an exponent h 
[Rypdal and Rypdal, 2010a, b]. Increasing h imphes higher 
persistence and predictability, corresponding to a reduction 
ol F. We compute h from the variogram (equation (12) in 
Zivkovic and Rypdal [2011]). Figure 3b shows no variation 
in ensemble mean of 2h for any of the indices, confirming 
the result obtained for F. 

The same analysis is done for B^, and v. Mean values for F 
and 2h are plotted in Figure 4, showing no significant change 
of predictabilty during substorms for the solar wind observ- 
ables. The example of two substorms from September, 8th, 
2002 shows a reduction of F in tAE around substorm onset 
(not shown here), but this day seems to be an exception 
rather than the rule. 

3.2. Determinism analysis 

The signals we study are dominated by a stochastic com- 
ponent and their L„ decreases when number of passes n 
is increased. However, the existence of a low-dimensional 
component in e.g. the tAE can be demonstrated by com- 
puting Ln and then do the same for the surrogate signals 
obtained by randomizing phases of Fourier coefficients. In 
order to compute mean L„ , we use AE index data from the 
year 2000; each L„ is computed over a record of 7200 points, 
such that the mean L„ is computed over N — 70 realizations. 
We have also computed mean L„ for several other lengths 
of records, but its value has not changed significantly. The 
embedding dimension is m = 6, and r = 20 minutes. The 
sample values of L„ are normally distributed, and hence the 
standard error of the ensemble mean is a/y/N. This error 
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Figure 7. L„ vs. n, a) squares are derived from tAL dur- 
ing substorms, triangles from the entire tAL time series, 
b) squares from tPC during substorms, triangles from 
entire tPC time series. Standard errors are of the same 
order as the symbol size. 



is for all L„-cmves shown in this paper of the same order as 
the symbol size, and is therefore not shown explicitly in the 
figures. As observed in Figure 5, I/„ for tAE is higher than 
for the randomized version, indicating existence of a low- 
dimensional and nonlinear component in the signal. The 
same is shown for tAL, tAU, and tPC indices. Notice that 
in case of fO-U in the inset of Figure 1, L„ for the random- 
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Figure 8. Ln vs. n, a) squares are derived from 
during substorms, triangles from entire B^ time series, 
b) squares are derived from v during substorms, trian- 
gles from entire v time series. Standard errors are of the 
same order as the symbol size. 
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Figure 9. a) Numerical solutions of equation (10), b) 
AE index. 
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ized version could not be distinguished from the Ln of the 
fO-U itself. 

Also, the same test was done for the geomagnetic index 
SYM-H in Zivkovic and Rypdal [2011], and even though 
analysis of this index yields low-dimensionality during mag- 
netic storms, L„ is indistinguishable from its randomized 
version when averaged over a year. Ln computed for 
and V is also indistinguishable from the L„ of their random- 
ized versions (not shown here). 

As we shall demonstrate below tAE, tAU and tPC all ex- 
hibit a discernible low-dimensional component during sub- 
storms. Since substorms occur very often, sometimes several 
per day, this low-dimensionality is discernible also when av- 
eraged over the year (as is shown in Figure 5). We use the 
substorm database and compute L„ in the time interval 1 
hour before and 1.5 hours after the substorm onset. This 
way we generally include all phases of the substorm. 

It has been shown in Zivkovic and Rypdal [2011] that 
the increase in the embedding dimension gives increased de- 
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Figure 10. Power spectral denities. Upper curve: from 
AE index. Middle: from solutions of equation (10) with 
w a fractional Gaussian noise with H = 0.45. Lower: 
from fO-U process with w the same as above. 
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Figure 11. Mean vector length Ln vs. number of pas- 
sages n: bullets are derived from low-pass filtered nu- 
merical solutions to equation (10), stars are derived from 
these solutions after randomization of phases of Fourier 
coefficients. Inset in the figure shows same characteris- 
tics for the filtered fO-U process. Standard errors are of 
the same order as the symbol size. 



terminism for the systems with low-dimensional dynamical 
component, while no change in determinism can be seen for 
the fO-U process. Therefore, the only constrain on the em- 
bedding dimension for our data should be the length of the 
studied time series. Ideally r should be chosen to corre- 
spond to the first minimum of the average mutual informa- 
tion (which is r « 20 in most of our data). However, in our 
data we are limited by the length of the time records, which 
is the typical length of a substorm. We use records of length 
2.5 hours (record length L = 150 data points), chose r = 10 
minutes, and m = 6 (when r = 20, results are the same, but 
due to better statistics, we use r = 10). This leaves a trajec- 
tory record in embedding space of length A*"— (m— 1)t = 100 
(from equation (2)), from which we compute a sample value 
of Ln for a given substorm. These values are computed over 
an ensemble of = 500 substorms; and the ensemble mean 
and standard error g are estimated. 

This analysis is made for tAE, tAL, tAU, and tPC, and 
compared to Ln computed from all the data from year 2000. 
Figures 6 and 7 illustrate that all geomagnetic indices (ex- 
cept tAL) exhibit elevated determinism during substorm 
times, i.e. low-dimensionality increases during substorms. 
Future studies should investigate why determinism in tAL 
index does not change during magnetospheric substorms. 
Notice that if we could obtain equally good statistics with 
higher embedding dimension, it would be very likely that the 
determinism during substorms would be even more elevated 
than it is for m = 6. 

Another way to demonstrate the existence of a low- 
dimensional component in these indices is to compare the 
determinism computed from the tAE during substorms 
with signals generated numerically from the fO-U equation, 
whose coefficients are fitted to the tAE under substorm 
through least square regression procedure. For each of a 
large number of realizations of the fO-U process, we compute 
the ensemble mean of L4 as a measure of low-dimensionality. 
The choice of L4 from the L„ vs. n curve is a compro- 
mise between clear separation between low-dimensional and 
stochastic dynamics and small error bars (which increase 
with increasing n). The ensemble mean for fO-U process 
is Li ~ 0.40 , while for deterministic process Li ~ 1. In 
contrast, Figure 6a shows that Z/4 ~ 0.65 for tAE under 
substorm conditions, which gives a clear indication that the 
AE index is neither fully deterministic nor stochastic. 

In Figure 8 we plot Ln for and v and observe no 
increase in determinsm during substorms. This indicates 




Figure 12. Mean vector length L„ vs. number of pas- 
sages n: bullets are derived from the entire low-pass fil- 
tered tAE time series, stars from this time series after 
randomization of phases of the Fourier coefficients. Stan- 
dard errors are of the same order as the symbol size. 
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that the increased low-dimensional component in the geo- 
magnetic indices during substorms is not imposed on the 
magnetosphere by the growth of such a component in the 
solar wind. In other words: the organization of the magne- 
tosphere during substorms is a seZforganization. 

4. Discussion 

We have applied recurrence-plot analysis and a test of 
determinism to geomagnetic indices AE, AU, AL and PC, 
as well as IMF and solar flow speed v. Recurrence 
plots were applied by March et al. [2005a] to connect solar 
wind observables to the AL and AU indices, concluding that 
the correlation between those indices results from magnetic 
storm signatures appearing in both time scries. Also, corre- 
lation between the solar wind electric field -iiBz and the AE 
index was studied in March et al. [2005b], invoking mutual 
information. Here it was found that correlation is present in- 
termittently on the timescales of a few hours, suggesting that 
substorms are information carriers between the solar wind 
and the magnetosphere. In our study, recurrence analysis is 
used to measure the inverse average diagonal line length T 
which is hcuristically shown to be a useful measure of pre- 
dictability of the dynamics. As an alternative predictability 
measure the self-similarity exponent h was measured during 
the course of substorms. No systematic variation of ft or F 
during substorms has been detected in an ensemble of five 
hundred substorms. This is true both for geomagnetic in- 
dices and the solar wind observables. Thus we conclude that 
geomagnetic indices which represent dynamics in the mag- 
netotail, plasma sheet boundary layer, partial ring current 
and polar cap convection, do not become more predictable 
during substorms. The same applies for observables repre- 
senting the dynamical state of the solar wind driver. From 
the analysis of a spatio-temporal dynamical model of the 
high latitude magnetic perturbation, it was shown in Val- 
divia et al. [1999] that the nonlinear dynamical model for the 
evolution of the spatial structure produces better prediction 
than the linear one during intense magnetospheric activity. 
In this article, a simple test for determinism has shown as 
well that the AE index exhibits some nonlinear characteris- 
tics, but in addition, it has s weak low-dimensional compo- 
nent. Those characteristics are elevated during substorms. 
Other geomagnetic indices, except the AL index, also show 
same signatures during magnetic substorms, which could in- 
dicate that the magnetosphere self-organizes and develops 
low-dimensional dynamics under substorm conditions. For 
the AE index similar results were obtained by Athanasiu 
and Pavlos [2001], who applied singular value decomposi- 
tion (SVD) analysis to the AE index. They concluded, by 
comparing the AE index with solutions of the Lorenz system 
contaminated by a colored noise term, that the first SVD 
component in both cases is entirely due to colored noise, 
while the higher-order SVD components are due to the inter- 
nal, low-dimensional and chaotic dynamics. By computing 
cross-correlation between the first SVD component of the 
AE index and the index itself, they concluded that the in- 
fluence of the flrst component is about 40 percent. Further, 
the first SVD component is attributed to the solar wind and 
is characterized as linear and stochastic, while higher SVD 
components arc attributed to magnetospheric dynamics. In 
our test of determinism the data arc not filtered to reduce 
the effect of the stochastic component on the analysis, but 
the results still reveal that unlike IMF Bz and solar wind 
flow speed v, the AE index contains components that makes 
it different from a high-dimensional or stochastic system. 

Thus, it seems firmly established that the auroral elec- 
trojet proxies exhibit an additional low-dimensional compo- 
nent associated with self-organization of the magnetosphere. 



the most prominent example being the auroral substorm. 
It should be stressed, however, that also for AE activity 
the major component is stochastic. Chapman and Watkins 
[2001] argue that the AE index contains information trans- 
ferred from the solar wind, since it has been shown in Free- 
man et al. [2000] that power laws in burst lifetime distri- 
bution for the AL and AU indices and for the solar wind 
electric field vBz are similar apart from a bump in the life- 
time distribution which Freeman et al. [2000] explain as a 
signature of the substorm current system. This indicates 
that the stochastic component of the auroral olectrojct ac- 
tivity to a great extent is a direct imprint of the solar wind 
turbulence. This conclusion is supported by the analysis in 
Rypdal and Rypdal [2010b], where Bz and tAE are found to 
exhibit very similar multifractal spectrum on the time scale 
up to 100 minutes. 

Our analysis supports a picture where the magnetosphere 
under quiet conditions resides in a forced state where the or- 
ganization of the magnetospheric dynamics reflects the or- 
ganization of the solar wind driver, i.e. the stochastic prop- 
erties of the global magnetospheric system and the driver 
are very similar. If substorms are trigged by solar wind fea- 
tures, this trigger is not an increase in organization or pre- 
dictability of the solar wind dynamics. During substorms ge- 
omagnetic indices are influenced by a self-organization that 
involves the major current systems in the magnotosphore- 
ionosphcrc, and hence indicates that a global instability is 
exited. This picture is consistent with the scenario described 
in Chang [1999], and conceptually simpler forerunners like 
Lewis [1991]. The latter sets out to explain the unpre- 
dictability of substorm onset, noting that external triggers, 
like reversal from northward to southward pointing Bz, are 
not always identifiable, and series of substorm cycles can oc- 
cur if southward IMF persists for prolonged times. South- 
ward IMF opens up for loading of magnetic flux and energy 
to the magnetosphere, and in Lewis [1991], and further elab- 
orated by Sitnov et al. [2001], parameters representing this 
loading are modeled as external time-dependent control pa- 
rameters, and the magnetospheric state-variable (e.g. a ge- 
omagnetic index) is modeled via a non-autonomous system 
on the form 

^=F{X,a,b), (9) 

where a{t), b{t) are two time-dependent external control pa- 
rameters. The underlying assumption is that the magne- 
tosphere resides in a forced equilibrium corresponding to a 
stable fixed point of this system, and hence these equilib- 
ria are located on the surface F{X, a,b) = in the three- 
dimensional (X, a, 6)-space. Choosing a third-order polyno- 
mial form, e.g. F(X, a, b) = -\- aX + b, gives rise to a 
folded surface that opens the possibility of a cusp catastro- 
phe which is interpreted as the substorm onset (expansion 
phase). According to this non-autonomous model for evo- 
lution of stable, forced equilibria the substorm onset is not 
really unpredictable since it depends on the control param- 
eters a and b. However, the catastrophe can occur along a 
curve in the (a, fe)-plane which may be hard to identify from 
this kind of conceptual model, and hence practical predic- 
tion may be difficult. 

One unsatisfactory feature of non-autonomous models of 
this kind is that the control parameters are not directly re- 
lated to the state of the solar wind driver, but rather a result 
of the solar- wind/magnetosphere interaction that is an in- 
tegral part of the dynamical system to be modeled. Thus, 
in this respect a more satisfactory approach is to model the 
magnetosphere-ionosphere as a dynamical system which is 
autonomous under constant forcing. Such models can be 
constructed with varying degree of sophistication, [Baker et 
al, 1990; Klimas et al., 1992; Horton and Doxas, 1998], and 
may give rise to chaotic signals that reproduce many of the 
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random characteristics of the magnctosphcric time scries. 
However, being deterministic and low-dimensional they can- 
not reproduce the strong stochastic component in the obser- 
vational signals. On the other hand, such models can be gen- 
eralized to include stochastic forcing from small-scale inter- 
nal dynamics and deterministic and stochastic forcing from 
variations of the driver. Thus, one may conceive complicated 
as well as simple conceptual dynamic-stochastic models that 
cam capture the essential stochastic dynamics as it presents 
itself in the observables studied here. In Rypdal and Rypdal 
[2010a, b] the AE index is modeled by a simple dynamic- 
stochastic equation that reproduces the general statistical 
features. The deterministic dynamics in this equation was 
represented by a drift term (a nonlinear damping) which 
prevents the solution to drift off to infinity, but the focus in 
those studies was on the time scales less than a 100 minutes, 
where the effect of the drift term is small. A version of this 
stochastic difference equation that exhibits on-off intermit- 
tency for certain choices of parameters is 

5X = M{x)5t + s/DX w (10) 

where M(x) = aX — X^ with a > models the drift term 
found from the AE index in Rypdal and Rypdal [2010b]. 
Here 5t is the discrete time step (the sampling interval of 
the time series) and SX{t) = X{t + St) — X{t). In Rypdal 
and Rypdal [2010a] w(t) is modeled as a particular multi- 
fractal stochastic noise process with unit variance. It was 
shown by Aumaitre [2005] that the on-off intcrmittcncy of 
this equation is sensitive to the nature of the noise term. 
If w{t) is approximated by a white Gaussian noise, equa- 
tion (10) in the limit (5—^0 reduces to the Ito stochastic 
differential equation 

dX = {aX - X^) dt + VD XdB{t), 

where B{t) is the Wiener process (Brownian motion). It 
can be shown from the associated Fokker-Planck equation 

that the stationary probability density for X is P{X) = 
(j^(2a/D)-i^-(x^/D)^ rpj^g divergence of P(X) as X -> for 
2o/-D < 1 is due to the on-off intcrmittcncy in this regime, 
which makes the solution reside in the vicinity of X = for 
a considerable portion of the time, while for 2a/ D— 1 small, 
but positive, the solution has an intermittent character more 
similar to the behavior of the AE index. Such a solution, 
for D = 0.1, is shown in Figure 9a, with a sample of the AE 
index shown in panel (b). Here w was chosen as a weakly 
anti-persistent fractional Gaussian noise with Hurst expo- 
nent H = 0.45, which is the H-value derived from the power 
spectral density for AE shown in Figure 10. The relation be- 
tween the spectral index /3, which is the slope of the straight 
lines fitted to the spectra plotted in a log-log plot, and the 
Hurst exponent of the differentiated signal is /3 = 2H + 1. 
On time scales < 100 time steps (minutes) AE index has 
/3 « 1.90, corresponding to 77 « 0.45. The power spectral 
density for the model signal shows a less clear power-law 
regime on these time scales, which is mainly due to a crude 
model for the nonlinear drift term M{X) in equation (10). A 
more "box-like" function would remedy this. On time scales 
> 100 mirmtcs the spectrum for the AE index has a pink- 
noise character (/? ~ 1), while the model time scries has a 
more gradual transition towards white noise. Further study 
on refinements of equation (10) is required to settle whether 
these spectral features are possible to reproduce within this 
class of one-dimensional stochastic equations. 

An interesting question is whether modeling along these 
lines can produce bursts (substorms) which enhance the de- 
terminism compared to a completely random process. As 
discussed before, the test based on comparing the I,„-curve 
with the one obtained after randomization of phases is a test 
on the existence of a low-dimensional and nonlinear compo- 
nent in the dynamics. Direct computation of L„ from the 



model signal docs not reveal a significant reduction of L„ af- 
ter randomization, in contrast to what was found from the 
MG system signal in Figure 1. One obvious reason for this is 
that the signal from equation (10) contains a strong stochas- 
tic component which is not present in the signal from the 
MG system. Hence, in this case it is necessary to perform 
a mild low-pass filtering of the model signal before com- 
putation of L„ [Kaplan and Glass, 1993]. The result after 
filtering is shown in Figure 11, indicating the presence of a 
low-dimensional nonlinear component. To make sure that 
this filtering does not introduce spurious low-dimensional 
nonlinearities in the signal, we perform the same test to a 
filtered fO-U process with similar spectral characterics as 
our model signal (the power spectral density for the fO-U 
signal is displayed in Figure 10). The filtered fO-U pro- 
cess shows very small change of Ln after randomization of 
phases, as shown in the inset of Figure 11. Obviously, low- 
pass filtering also have an effect on this test when applied to 
physical signals with a stochastic component. In Figure 12 
we show the effect of applying the same filter to tAE, which 
should be compared to the result for the unfiltered signal 
shown in Figure 5. 

The nonlinearity producing determinism in the model sig- 
nal from equation (10) is a combination of the nonlinear 
drift term and the multiplicative noise term -s/DXw. The 
multiplicative term can be eliminated by the transformation 
Y = logX, but then Ito's formula [Gardiner, 1985] yields a 
new drift term on the form M{Y) = e'^ M{e^) - D/2. 
Note that if our original drift term is a linear damping 
M{X) = —vX the transformed drift term reduces to a nega- 
tive constant M{Y) = -{^+0/2). This yields Y{t) -oo, 
and hence X{t) Q as t ^ oo, and demonstrates that the 
drift term must be nonlinear to produce stationary time se- 
ries from a model with a multiplicative noise term like equa- 
tion (10). 
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